Tutorial Overview
Biodiversity is a hot topic in the fact of unprecedented global change in the 21st century. With a publicly available dataset, from the Turkey Lakes Watershed, we demonstrate how to quantify biodiversity quantitatively and visually using R.
Learning objectives
Utilizing historical data from a freshwater lake ecosystem, by the end of this tutorial you will be able to:
- Identify, calculate, and implement commonly used metrics to quantify biodiversity
- Understand and create figures uses to display biodiversity data
- Apply this concepts to other aspects of the dataset used in this tutorial
Ecologists are often interested in the effects of perturbations (e.g. fires, logging, flooding or pollution) on biodiversity, which can yield information and context for ecosystem health. However, boiling biodiversity down to a single metric can be difficult and incomplete depending on the goal. With the help of the Turkey lake dataset of benthic invertebrates, we will present you with a number of common metrics to assess biodiversity, their meaning, and how they can be visualized with the help of Rank-Abundance curves.
You will then have the option to learn to calculate biodiversity metrics yourself using an R package called codyn.
Why biodiversity matters
In an ever-changing world, we are often concerned with how organisms will be able to cope with disturbances and increasingly, global climate change. Biodiversity is the variety of life at all levels, from genes to ecosystems. Not only is biodiversity important for ecosystems to function properly, it is also an important buffer to withstanding changing environments. When there is high biodiversity, there is a higher chance of some organisms surviving ecosystems that are changing due to disturbance or climate change. Biodiversity therefore increases the resiliency of these ecosystems and the likelihood that they will be able to persist beyond this change.
Types of data typically collected
Biodiversity is a broad term, and can be difficult to quantify. Biologists have created several standard measures of biodiversity that convey different information about ecosystems and the species living there.
Often times, biologists interested in collecting biodiversity data collect a standard set of responses. These often include:
- The different kinds of species in a given area or environment
- The number of individuals of each species (abundance)
- Information about the area itself (i.e. habitat type, weather conditions, etc.)
- If the work involves an setting up an experiment, collecting the control vs. experimental groups as well as metrics for replication and statistical power
- The time frame over which the data were collected
Ways to calculate diversity
Once we have collected our biodiversity data, there are several different ways to analyse and understand changes in community composition, or the relative abundances of all taxa in a community. Each of these metrics allow us to look at different aspects of the biodiversity data we collect.
- Species richness quantifies the number of different species in a community.
- The Shannon-Weiner Index estimates species diversity by taking into account the number of species living in a habitat (richness) and their relative abundance (evenness).
- Simpson’s Diversity Index also takes into account species richness and evenness, but puts more emphasis on rare species.
- The Evenness Index assesses the evenness of the abundances of each species.
The Data
To explore these different diversity metrics and how they are calculated, we have simulated a biodiversity dataset with 10 unique species and random abundances across 3 years. This will allow us to look at how diversity changes over time.
Here, we generated a simulated dataset of different abundances of species across 3 years. Each number in the species columnn represents a unique species, abundance is the number of individuals counted in a given year. Here is the first 6 rows in our simluated dataset.
hide
pacman::p_load(tidyverse, patchwork)
set.seed(222)
sim1 <- data.frame(
species = sample(1:10, 20, replace=2),
abundance = sample(0:200, 20, replace=2),
year= 2001)
set.seed(723)
sim2 <- data.frame(
species = sample(1:10, 20, replace=2),
abundance = sample(0:200, 20, replace=2),
year= 2009)
set.seed(117)
sim3 <- data.frame(
species = sample(1:7, 20, replace=2),
abundance = sample(0:200, 20, replace=2),
year= 2020)
sim <- rbind(sim1, sim2, sim3)
head(sim) species abundance year
1 2 149 2001
2 7 104 2001
3 8 42 2001
4 6 8 2001
5 10 191 2001
6 4 196 2001Here is a look at the simulated data.
hide
hide
#facet_wrap(~year)Figure 1. Simulated species abundances for 3 different years and 10 different species. Bar are the abundances of the species, coloured by year.
We can see that this data is quite messy and hard to find clear patterns in changes across species, abundance, and years. Below are metrics we can use to better define and clarify these patterns.
Species richness
Perhaps the simplest measure of biodiversity is species richness (\(S\)).
The equation for species richness is:
\(S = \sum_{i=1}^{S}p_i^o\)
(Describe what the variables mean).
Step 1: Count the number of unique species recorded in that year. Then you can look at how the number of different species in a community changes over time.
Interpretation: A higher number means there are more species, and therefore species richness is higher.
hide
Error in `group_by()`:
! Must group by variables found in `.data`.
âś– Column `time` is not found.hide
#we should probably make this into a nice looking kable?
kable(richness_table, col.names=c("Year", "Richness"), align="c")%>%
kable_classic(full_width = F, html_font = "Cambria")Error in kable(richness_table, col.names = c("Year", "Richness"), align = "c"): object 'richness_table' not foundWe can see in this table that across years, the number of species in the same area fluctuated, some years with more species, other years with less.
Shannon-Wiener Index (H)
Another metric we can use to analyze biodiversity is the Shannon-Wiener Index (\(H\)), which indexes across the sum of all relative proportions of species richness’s with an additional logarithmic transformation. This equation uses both the number of species in a specific area (richness) and relative abundance (evenness).
The equation for the Shannon-Weiner Index is:
\(H = - \sum_{i=1}^{S}p_ilnp_i\)
Step 1: Calculate \(p_i\), which is the proportion of \(i\)th individuals in the community. \(p_i = n_i/N\)
Where \(n_i\) is the number of individuals in the ith species and and \(N\) is the total number of individuals in the community.
Step 2: Multiply the proportion of each species (\(p_i\)) by the logarithm of the proportion (\(lnp_i\))
Step 3: Sum each of these values for each species.
Step 5: Multiply the sum by -1.
Interpretation: The minimum value of the Shannon’s diversity index is 0, which means that there is no diversity (i.e. only one species is found in that habitat). The values increase as the number of species increase, and is maximized at a given number of species when is an equal abundance of each species.
hide
Error in `group_by()`:
! Must group by variables found in `.data`.
âś– Column `time` is not found.hide
kable(shannon, col.names=c("Year", "Shannon's Index"), align="c")%>%
kable_classic(full_width = F, html_font = "Cambria")Error in kable(shannon, col.names = c("Year", "Shannon's Index"), align = "c"): object 'shannon' not foundFrom our simulated dataset, we can see that the Shannon diversity increases over time.
Simpson’s Index (D)
Another metric we can use to analyse biodiversity in the Simpson’s index (\(D\)) which indexes across the sum of all relative proportions of species’s richness with an additional square power transformation. Though very similar to the calculation in the Shannon-Wiener Index, the Simpson index is more focused on dominance of species as it highlights the proportion of species in a sample.
The equation for Simpson’s diversity is:
\(D = 1-(\sum n_i(n_i-1))/N(N-1))\)
Where \(n\) is the number of individuals of the ith species and \(N\) is the total number of individuals in the community.
Step 1: Multiply the number of individuals of a given species (\(n_i\)) by (\(n_i-1\)).
Step 2 Multiply the total number of individuals in the community (\(N\)) by (\(N-1\)).
Step 3: Divide the number from step 1 by the number from step 2.
Step 4: Once you have a number from step 3 for each species in the community, sum all of these together.
Interpretation: The higher the value of \(D\), the greater the diversity in the community. An index close to 1 means that there are several species in the community, and the population proportion of species is even.
hide
Error in `group_by()`:
! Must group by variables found in `.data`.
âś– Column `time` is not found.hide
kable(simpson, col.names=c("Year", "Simpson's Index"), align="c")%>%
kable_classic(full_width = F, html_font = "Cambria")Error in kable(simpson, col.names = c("Year", "Simpson's Index"), align = "c"): object 'simpson' not foundThe Simpson index increases over time, which tells us that by this measure, biodiversity is increasing.
Evenness Index (E)
Additionally, we can think about the evenness of species across a given area. Evenness is a metric for assessing species dominance. If evenness is high, it means most species are of equal dominance. If evenness if low, it means some species are more dominant (i.e. have higher abundance) than others. Evenness is maximized when all species have equal abundances.
The equation for species evenness is:
\(E = H / ln(S)\)
Step 1: Divide the value of the Shannon-Wiener index (\(H\)) by the species richness (\(S\)).
Interpretation: Values closer to one signify more evenness and values closer to 0 signify lower evenness.
hide
Error in `group_by()`:
! Must group by variables found in `.data`.
âś– Column `time` is not found.hide
kable(evenness, col.names=c("Year", "Evenness"), align="c")%>%
kable_classic(full_width = F, html_font = "Cambria")Error in kable(evenness, col.names = c("Year", "Evenness"), align = "c"): object 'evenness' not foundIn our simulated dataset, we can see that evenness is low from the start and increases over time, suggesting that species abundances are becoming more even.
Visualizing species diversity
Rank abundance curves with simulated data
One way to visualize biodiversity data is to make a Rank Abundance Curve (RAC). On the y-axis we have the species (numbered), and on the y axis we have the abundance. The plots are organized from high to low abundance. Using RACs, we can visualize the change in relative abundance and species dominance over time.
hide
library(tidyverse, viridis)
# Code or exercise on calculating species richness?
#before
pacman::p_load(tidyverse, patchwork)
data_2007 <- data.frame(
species = c(1,2,3,4,5,6),
abundance = c(77, 109, 24, 38, 11, 6),
time= 2007)
data_2018 <- data.frame(
species = c(1,2,3,4,5,6),
abundance = c(47, 79, 2, 15, 169, 111),
time= 2018)
data <- rbind(data_2007, data_2018) %>%
group_by(species, time) %>%
summarise(abundance= sum(abundance))
#2007
a<- data %>%
filter(time==2007) %>%
ggplot( aes(x=reorder(species, -abundance), y=abundance, shape=as.factor(species), colour=as.factor(species)))+
geom_point(size=4)+
scale_colour_viridis_d()+
ggtitle("2007")+
labs(x="Species", y="Abundance")+
theme_classic()+
theme(legend.position="none")+
labs(x="Species", y="Abundance", color="Species")
#2007
b<- data %>%
filter(time==2018) %>%
ggplot(aes(x=reorder(species, -abundance), y=abundance, shape=as.factor(species), colour=as.factor(species)))+
geom_point(size=4)+
scale_colour_viridis_d()+
ggtitle("2018")+
labs(x="Species", y="Abundance")+
theme_classic()+
theme(legend.position="none")+
labs(x="Species", y="Abundance", color="Species")
a / b
We can see here between 2007 and 2018 in our simulated community, there were differences in absolute abundance as well as which species were dominant, which species disappeared, and which species were new. For instance, in 2007, species 2 was the most abundant species, but it 2018, it was the 3rd most abundant species. However, species 3 is lost by 2018, while species 6 is gained in 2018.
This where I couldn’t seem to get the second panel set to show up
Rank abundance curves with Turkey lakes invertebrate data
The Data
The Turkey Lakes watershed is located in Ontario, Canada, approximately 50km North of Sault Ste. Marie. It is 10\(km^2\) in area and contains a chain of 4 lakes. The Turkey Lakes Watershed Study (see map) is an initiative by several agencies of the government of Canada, initially designed to study the effects of acid rain on forest ecosystems. From 1995 to 2009, scientists collected, identified and counted benthic invertebrates from various stream beds around the Turkey Lakes catchment.
Benthic invertebrates are small, often microscopic organisms (see Fig. 1), but form an important link between aquatic and terrestrial habitats . They can decompose terrestrial matter such as leaves, or consume periphyton (algae and cyanobacteria) growing on rocks within their streams. Benthic invertebrates can be a food source for aquatic animals, like small fish, or terrestrial animals such as birds. An experiment was conducted in 1997 where certain sampling sites experienced different levels of logging intensity. Due to their interaction with the terrestrial realm, we can hypothesize that the biodiversity of benthic invertebrate communities might be affected by logging.
Figure 2. Interactive map of the Turkey Lakes study site
Load Data
These data are from the Turkey Lakes watershed experiment. This dataset contains the abundances of benthic invertebrate species measured in the May and June from 1995 to 2001. They recorded the type and number of invertebrates in an area of the stream, with 10 replicates per sampling event. These forests next to these streams were logged in 1997 with varying harvest intensities, with low, medium, and high intensity harvest and a control stream with no forest harvest. We will using the codyn package to analyze the diversity of these benthic invertebrates over time.
year species treatment replicate abundance
1 1995 Alloperla control 1 0
2 1995 Alloperla control 2 0
3 1995 Alloperla control 3 4
4 1995 Alloperla control 4 3
5 1995 Alloperla control 5 3
6 1995 Alloperla control 6 3Raw data
RACs in nature are rarely as evenly distributed as the simulated data. There is often few highly abundant species and many very rare species. Below, you can see the rank-abundance distribution from one of the catchments in our Turkey Lakes benthic invertebrate dataset from the year following logging treatments (1998).
Log transformed data
To better visualize the distribution, the y-axis is typically log-transformed
Rank
Because most classic metrics of biodiversity do not take the species identity into account, we can drop species names from the x-axis, and simply refer to it’s rank in the community. This is where the term rank-abundance comes from.
Curve
Finally, instead of using bars to represent abundances, we can simply use lines to trace the distribution. This allows us to compare multiple distributions, in this case, from different three different sites with three levels of logging intensity.
Species Richness
We can compare species richness by comparing the maximum rank.
hide
dfRank = df %>%
filter(TotalCount > 0)
annotations = dfRank %>% summarise(maxRank = max(Rank))
annotations$xend = c(40, 31, 34)
annotations$yend = rep(1, 3)
annotations$x = annotations$xend + 3
annotations$y = rep(5, 3)
ComparisonPlot +
geom_segment(data = annotations, aes(yend=yend, xend=xend, x=x, y=y), size=1.2, arrow = arrow(), show.legend = FALSE) +
geom_text(data = annotations, aes(x = x+1, y=y+2, label=maxRank), size=6, show.legend = FALSE) 
Eveness
The biodiversity metrics described (Shannon, Simpson, etc.) also incorporate community evenness, which can be thought of as the slope of rank-abundance curve. For example, a community with one very dominant species and a rapid drop-off in species abundance will have a lower evenness than a community where more species are common. In the plot below, the community at the site of “medium” logging (orange) is less even than the other two since it has a steeper average slope.
hide
q = {ComparisonPlot + geom_smooth(method=lm, se=FALSE, data=dfRank) }%>%
ggplot_build()
q$data[[1]]$alpha = 0.5
ggplot_gtable(q) %>%
plot 
Several models have been developed to understand rank-abundance curves (e.g. log-normal, broken stick model).
The radfit function from the vegan R package allows fitting of these different curve models (see https://search.r-project.org/CRAN/refmans/vegan/html/radfit.html)
codyn tutorial
Load Packages
We will be using the codyn package for this tutorial. We will also use ggplot2 and dplyr for plotting and cleaning.
For many functions, the codyn package can only look at data from one stream at a time. So first we will separate each stream into it’s own data frame.
hide
# high intensity logging treatment
subset_codyn_high <- data %>%
filter(treatment == "high")
# medium intensity logging treatment
subset_codyn_med <- data %>%
filter(treatment == "medium")
# low intensity logging treatment
subset_codyn_low <- data %>%
filter(treatment == "low")
# control stream (no logging)
subset_codyn_ctrl <- data %>%
filter(treatment=="control")Using the codyn package to analyze diversity data over time
Ecologists are often interested in how diversity changes over time in one place. This is particularly interesting when you are looking at how species diversity responds to environmental change. This environmental change can be a single disturbance event such as a forest fire or more long-term change such as temperature increases due to climate change. The Turkey Lakes dataset provides us with an opportunity to look at how the diversity of benthic invertebrates changes over time and in response to logging. We can compare the diversity of different streams with different logging intensities and we can also compare the diversity of streams.
The codyn package allows us to analyze the change in biodiversity over time using different metrics.
Using community_structure()
community_structure() is a function in the package codyn that calculates species richness and evenness.
Species Richness and Evenness
For example, we can use community_structure() to look at the change in species richness and evenness over time in one of the Turkey Lake streams. We will choose the stream with high logging intensity as an example. This function can only look at one stream at a time so make sure to subset out the stream that you are interested in.
Let’s take a look at the high intensity logging stream data.
hide
head(subset_codyn_high) year species treatment replicate abundance
1 1995 Amphinemura high 1 9
2 1995 Amphinemura high 2 24
3 1995 Amphinemura high 3 53
4 1995 Amphinemura high 4 15
5 1995 Amphinemura high 5 59
6 1995 Amphinemura high 6 37As you can see, there is a column for year, species, replicate, and abundance.
The community structure function has five arguments.
1. df is where you put in the data frame of interest, in our case subset_codyn_high.
2. The time.var argument is where you put in your time variable (year)
3. The abundance.var argument is where you put the column that holds the abundance variable (abundance).
4. The replicate.var argument holds the column with the replicates (replicate)
5. The metric argument allows you to specify the metric of evenness that you want the function to calculate. There are three possible metrics: Evar, SimpsonEvenness, and EQ, which all calculate evenness in different ways. We will use Simpson’s Evenness in this tutorial.
Remember, if you ever need a refresher on the syntax you can type ?community_structure into the R script.
hide
richness_high<-community_structure(df=subset_codyn_high, time.var="year", abundance.var="abundance", replicate.var="replicate", metric="SimpsonEvenness")We can take a look at the output.
hide
head(richness_high) year replicate richness SimpsonEvenness
1 1995 1 17 0.2210248
2 1996 1 11 0.5785671
3 1997 1 16 0.3436584
4 1998 1 18 0.2786702
5 1999 1 15 0.3540228
6 2001 1 23 0.2661854As you can see, this returns a data frame with the year, replicate, richness calculation, and evenness measure. It calculates richness and evenness separately for each replicate sample in the stream, so we have a measure of variability within samples as well as variability across years.
Plot species richness
We can plot this data frame to see the changes over time.
hide
ggplot()+
geom_point(richness_high, mapping=aes(x=year, y=richness))+
theme_classic()+
labs(y="Species Richness", x="Year")+
geom_vline(xintercept=1997, linetype=2, color="red")
Figure Species richness of aquatic macroinvertebrates over time. Each point is a replicate.
Here is a plot of species richness, with a different value for each replicate. The red line is the year that the streams were logged. Do you see any change? Do you think there are any significant patterns? How does the variability between replicates compare to the variability between years?
Plot species evenness
We can also plot species evenness over time.
hide
ggplot()+
geom_point(richness_high, mapping=aes(x=year, y=SimpsonEvenness))+
theme_classic()+
labs(y="Simpson's Evenness", x="Year")+
geom_vline(xintercept=1997, linetype=2, color="red")
Figure Species evenness of aquatic macroinvertebrates over time. Each point is a replicate.
Do you see any patterns in the evenness between years?
Using community_diversity()
Shannon’s diversity
Species richness and evenness are very simple ways to look at species diversity. However, as you learned from the rank abundance curves, species have different abundances within the community, with some species having very high abundance and most species having low abundance. Species richness does not take into account the relative abundances of each species, which limits the inferences we can make. The Shannon-Weiner diversity metric does take into account these relative abundances when calculating the diversity of a community and allows us to make more meaningful inferences about how biodiversity is changing over time.
The codyn package allows you to calculate the Shannon diversity index for a single site over time using the community_diversity() function. We will again use the data from the stream with the high logging intensity. Luckily, this function has the exact same syntax as the community_structure() function! The only difference is that the options for the metric argument are Shannon and InverseSimpson. We will first use Shannon.
hide
shannon_div_high<-community_diversity(df=subset_codyn_high, time.var="year", abundance.var="abundance", replicate.var="replicate", metric="Shannon")Let’s take a look at the output for this function.
hide
head(shannon_div_high) year replicate Shannon
1 1995 1 1.836002
2 1996 1 2.044210
3 1997 1 2.107977
4 1998 1 2.016146
5 1999 1 1.948510
6 2001 1 2.297970You can see it gives us a data frame of year, replicate, and Shannon, which is our measure of Shannon diversity.
Plot Shannon’s diversity
We can plot these to see change over time.
hide
ggplot()+
geom_point(shannon_div_high, mapping=aes(x=year, y=Shannon))+
theme_classic()+
labs(y="Shannon Diversity", x="Year")+
geom_vline(xintercept=1997, linetype=2, color="red")
Figure Shannon’s diversity of aquatic macroinvertebrates over time. Each point is a replicate.
Simpson’s diversity
We can also use the same function to calculate Simpson’s diversity by simply changing the metric argument to InverseSimpson.
hide
simpson_div_high<-community_diversity(df=subset_codyn_high, time.var="year", abundance.var="abundance", replicate.var="replicate", metric="InverseSimpson")We can look at the output.
hide
head(simpson_div_high) year replicate InverseSimpson
1 1995 1 3.757421
2 1996 1 6.364238
3 1997 1 5.498534
4 1998 1 5.016063
5 1999 1 5.310341
6 2001 1 6.122265Plot Simpson’s diversity
And then plot.
hide
ggplot()+
geom_point(simpson_div_high, mapping=aes(x=year, y=InverseSimpson))+
theme_classic()+
labs(y="Simpson Diversity", x="Year")+
geom_vline(xintercept=1997, linetype=2, color="red")
Figure Simpson’s diversity of aquatic macroinvertebrates over time. Each point is a replicate.
How do these patterns compare to the patterns of Shannon’s diversity over time? How do the values compare?
Activity
Calculate richness, evenness, and diversity using both the Shannon and Simpson indices for the streams with low intensity logging, high intensity logging, and no logging. The data for these streams are subset_codyn_low, subset_codyn_high, and subset_codyn_ctrl. Plot these indices. Do you see any patterns across years? How do these compare to the patterns in the high intensity logging stream?
References
I got the equations from this paper: DOI:10.1007/s10531-016-1261-0